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ABSTRACT - A decentralized control framework is investigated for applicability 
of formation flying control in libration orbits. The decentralized approach, being 
non-hierarchical, processes only direct measurement data, in parallel with the 
other spacecraft. Control is accomplished via linearization about a reference 
libration orbit with standard control using a Linear Quadratic Regulator (LQR) 
or the GSFC control algorithm. Both are linearized about the current state 
estimate as with the extended Kalman filter. Based on this preliminary work, the 
decentralized approach appears to be feasible for upcoming libration missions 
using distributed spacecraft. 

1 - INTRODUCTION AND BACKGROUND 

Libration orbit missions in the near future will use multiple spacecraft in a distributed approach to 
perform interferometry and optical measurements not achievable by single spacecraft. We 
investigate formation flying concepts in a Sun-Earth libration orbit dynamical environment for 
these future missions. The Guidance, Navigation, and Control Center (GNCC) at the Goddard 
Space Flight Center (GSFC) is developing and implementing enhanced autonomous formation 
flying systems which improve operations while minimizing impacts to onboard software and 
hardware requirements, software development, and integration ([Folt98], [Carp99], [CarpOO]). In 
future missions, distributed systems are expected to require a communications cross-link for sharing 
housekeeping telemetry, navigation data, and measurements. In anticipation of such missions, the 
GSFC GNCC is investigating an option for closed-loop autonomous navigation and maneuver 
control of satellite formations that are based on the decentralized framework developed in [Spey79]. 
We begin by providing a brief description of libration orbit simulation, then describe GSFC 
decentralized control methods useful for support of these missions. We demonstrate new 
developments such as GSFC’s formation control algorithm and traditional linear quadratic control 
methods which are incorporated into the decentralized approach. 

Problems of single spacecraft control of circular restricted three-body (CRTB) motion have been 
previously investigated using state-space equations to characterize the linearized equations of 
motion ([Hoff93], [Farq70], [Wie98]). State-space analysis methods in control theory provide a 
useful framework for defining goals and the optimal control of satellites designed to fly about a 
reference orbit. We advance formation flying control using state-space by incorporating two control 
algorithms into a decentralized architecture. The mathematical foundation for the explanation of 
CRTB motion is briefly addressed, but only to the degree required to understand the results of the 
simulations presented. We consider the effects of third body and solar radiation perturbations into 
the state space simulations and calculate command responses, including disturbance rejection, on 
the formation of two spacecraft. 

The decentralized approach is non-hierarchical, eliminating the need tor coordination by a central 
supervisor and permits graceful degradation of system performance in the presence of detected 




failures. Each spacecraft (spacecraft) in the decentralized network processes only its own 
measurement data, in parallel with the other spacecraft. Although the total computational burden 
over the entire network is greater than it would be tor a single, centralized controller, fewer 
computations are required locally at any individual spacecraft. Requirements tor data transmission 
between spacecraft are minimized, at the cost ot locally maintaining an additional data vector. This 
data vector retains a memory of all past measurement history from all of the spacecraft compressed 
into a single vector with the dimension of the state. 

2 - LIBRATION ORBITS 

A wealth of references can be readily found on methods to derive the principle equation for the 
dynamic behavior of a single spacecraft in orbit around the Sun-Earth co-linear libration point, L 1 . 
These methods using Hamiltonian, Jacobian or Newtonian approaches start with the same initial 
assumptions: an infinitely small mass is placed in the gravitational field of two massive bodies or 
primary bodies. These methods are generally referred to as the restricted three-body (RTB) 
problem. Furthermore, if it is assumed that the primary bodies follows a circular orbit around their 
barycenter, then this is called the circular restricted three-body (CRTB) problem. An appropriate 
coordinate system needs to be selected. 

2.1 - A State Space Model 

As with the current literature about this subject ([Hoff93] [Wie98]), we use three-dimensional 
Cartesian coordinates denoted by capital letters to describe a system with an origin at the barycenter 
of the primary bodies and small letters will be used for an origin at a libration point. Figure 1 
illustrates the coordinate system and sample orbit that will be used here. Here m is an infinitely 
small mass in the gravitational field of the primary bodies M\ and M 2 . 



Fig. 1: Coordinate system and sample libration orbit 


The linearized equation of motion for m close to the libration point is 

x-2ny = U xx x , y + 2/rc = t/ n ._y , and z=U zz z (2.1) 

where U rY =—, U YY and U u =^. Uxx, U YY and U zz are calculated at the respective 

.tv dx2 YY qy- dZ- 

libration point to get the respective equation of motion and are constants. As a result of the 
linearization, x and y are coupled whereas z is now completely independent and is a simple 
harmonic. These equations can also be written in state-space form as below where x / represents the 
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state of the yth spacecraft in the formation. This is the torm that will be used in the following 
sections. 

2.3 - Reference orbit 

In order to compare our simulations using the state space matrix to a reterence and to have a 
reference orbit to control a spacecraft to, it is essential to find a proper reference model. The model 
used generally in the literature has an oscillating form. The amplitudes and phases of this reterence 
model depend on the initial conditions of the simulation. 

x = - A x sin (<o xy t + <p xy ) x = - A x co xy cos(co xy t + <p xy ) 
y = -A y cos((d xy t+(p xy ) y = A y co xy sin(co xv t + (p xy ) (2-3) 

z = A, sin (a)J+<p. ) i = A.co. cos (coj +<p,) 

2.4 - Simulation: The Circular Restricted Three-Body Problem for Earth-Sun LI 

To visualize the results above in graphs, a simulation is implemented with MATLAB/SIMULINK 
that will calculate the motion and plot the results. To achieve more accurate results in a CRTB, 
the Earth system will consist of the mass of Earth plus the mass of moon. This assumption is 
generally accepted by [Sze67], [Hoff93] and [Wie98]. Using the equations above it can the be 
computed that: 

M l =1.9891xl0 30 kg, Mi =6.0477x1 0 24 kg (Mass of Sun and Moon) 

Uxx= 3. 6262x1 O' 13 m/s 2 , U n =- 1.2 185x1 O' 13 m/s 2 ' f/zz=- 1.6 149x1 O' 13 m/s 2 

coxy = 4.1 603x 1 O' 7 rad/s* Txy= 174.8 days and co z = 4.01 85 5x1 O' 7 rad/s » T, = 180 days 

With these values a libration orbit is propagated over five years. Figure 2 shows the result of this 
simulation. The simulation error compared to the reference model if expressed as a percentage of 
the absolute distance in each direction, would not exceed 0.067%. The solar radiation pressure and 
the effects from the Moon and Jupiter were not included in the formation control aspect of this 
paper, but should be considered in future analysis. 



Fie. 1: Libration orbit around Sun-Earth LI without nerturbations 

3 - CONTROL ALGORITHMS 

The analysis presented here is based on generalizing earlier work [Hoff93] [Wie98] to the task of 
distributed spacecraft at the Sun-Earth, rather than Earth-Moon, collinear equilibrium point. In 
[CarpOO], the decentralized LQG framework developed [Spey79] is applied to an idealized (linear, 







time-invariant) version of a low Earth orbit satellite formation control problem. In this paper, the 
approach taken in [CarpOO] is applied to a centralized design for a formation of satellites at the Sun- 
Earth L\ point, which is based on the work of [Hoff93] and [Wie98]. Also, the Folta-Quinn (FQ) 
orbit control algorithm of [Folt98] is substituted for the LQR control used by [CarpOO], as was done 
in [Carp99], This is accomplished by treating the framework of [Spey79] as a decentralized 
estimator only; globally optimal states, rather than globally optimal controls, are then reconstructed 
at each spacecraft for use by the FQ algorithm. 

3.1 — Centralized LQR Design 

To maintain a spacecraft in orbit about L\ a control can be introduced that provides the propulsion 
system the direction and magnitude of the thrust in order for it to stay as close as possible to a 
reference trajectory. One such control used by spacecraft is the Linear Quadratic Regulator (LQR). 
The LQR is a control algorithm that uses state feedback to keep a system like a spacecraft to an 
equilibrium condition and then holds it there despite disturbances. State feedback indicates that the 
control command is directly proportional to the state. In the case of a spacecraft, the control 
command would be the thrust and the state would be the position and velocity of the spacecraft 
relative to the libration point. The reference orbit would define the equilibrium condition [Wie98]. 

The state dynamics of the controlled spacecraft are described by 

x' =A J x J +B'u' +»,„i pn ,, + a (3-1) 

where the matrix & maps the control input from the control space to the state space. To keep the 
spacecraft as close to the reference orbit, x R , as possible, and to minimize fuel usage, the 
performance index 

J =— J{(x-x") f g(x-x /f )+u / /?u}c?/ (3.2) 

is minimized. To control a formation of satellites, two approaches are possible. Either each vehicle 
can be individually controlled to its own reference trajectory about LI, or the spacecraft relative 
positions can be controlled. In the former case, the control problem can be decoupled such that 
each vehicle’s optimal control is independent of the other spacecraft states and controls. Then, x 
and u above represent a single spacecraft state and control x 7 and u 7 . Otherwise, x and u above 
represent all the spacecraft states and controls, x = [x 1 ; x 2 ; ...; x A ], u = [u 1 ; u 2 ; ...; u ]. In either 
case, the solution to this problem is 

u'=-M~V] 7 Sx (3.3) 

where S is the steady-state solution of the algebraic Riccati equation for time invariant problems 

S{BR-‘B' )S-SA -A'S-Q=0 

The solution to this equation is calculated using the MATLAB/Control Systems Toolbox function 
Iqr. The Iqr function needs the A, B, Q, and R matrices as its input. The matrix A has been defined 
previously for a single spacecraft. For multiple spacecraft, A is assumed to be equivalent for each 
spacecraft as the separation distance between spacecraft for our analysis is much smaller that the 
distance from either the Earth or Sun. If all the spacecraft are to be simultaneously controlled, A is 
defined to include a copy of the single spacecraft A 7 matrix for each vehicle in the formation along 
its diagonal. 

The B matrix maps the control input from the control space to the state space, e.g. B determines 
which states are influenced directly by the control. For the centralized controller design model it 
will be assumed that the spacecraft thrusts continuously and that the thrust level can vary 
continuously in time. Also, the spacecraft’s control system can apply the thrust directly in each 
component of the coordinate system [Wie98]. Therefore, the B matrix for an individual spacecraft is 
given by 



y 



(3.4) 



If all the spacecraft are to be simultaneously controlled, B is defined to include a copy ot the single 
vehicle's B matrix for each vehicle in the formation along its diagonal. 

The O and R matrix are performance weights for the LQR. O is the weighting matrix tor the state 
error7which determines how accurate the reference is followed, and R is the weighting matrix tor 
the control effort, which determines the damping of the system. In both cases trial and error and a 
certain amount of experience determine the selection of O and R. There are however some 
guidelines proposed by [Wie98] on how to proceed. Let 



The value of r determines how well the controlled system is damped. Ifp and q are selected, then r 
can be used to tune the controller until it produces reasonable results. In the present case, p = 100 
m 2 , corresponding to a position tolerance of 10 m, and q = 0.0001 (m/s) corresponding to a 
velocity tolerance of 0.01 m/s. If all the spacecraft are to be simultaneously controlled, Q and R are 
defined to include copies of the single vehicle’s Q and # matrices for each vehicle in the formation 
along their diagonals. 

3.2 - Disturbance Accommodation Model 

The spacecraft thrusts continuously due to the perturbing forces and because the reference model 
does not correspond with the numeric solution exactly. To save fuel it is however desirable to 
decrease thrusting time. [Hoff93] and [Wie98] have similar methods to limit fuel consumption by 
introducing a disturbance model. This technique was successfully used to improve the attitude 
control design for the International Space Station. 

The reason that the LQR design cannot keep the thrust at a lower level is that the A matrix does not 
include the disturbances from the perturbations and the fact that the simulation does not equal the 
reference model exactly. The real solution has additional nonlinear terms that cause the controller to 
signal a higher thrust. A disturbance model can reduce the effects of these terms and reduce the 
thrust. The disturbance accommodation model allows the states to have non-zero variations from 
the reference in response to the perturbations without inducing additional control effort. In essence, 
the additional states introduced by the disturbance model are in resonance with modes of the control 
effort caused by the disturbances. As the disturbance states are stimulated, they absorb control 
effort, so that in steady-state, the control effort goes to zero. 

The state matrix for the disturbance model can be written as 

8 7 =4LA' ( 3 - 6) 
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where o' - [a, and & = [0, represent n periodic disturbances with frequencies given 
bv to 7 = ;cd„], and x 7 = [T,;...;x m ] represent m constant disturbances. To completely 

accommodate all disturbances the model would have to include many periodic disturbances and 
integral states. This however would make the model very complex and computational expensive. 
Therefore, only the most severe disturbances are included in this model. 



The task at hand is to find a suitable set of frequencies. A method that is described in [Wie98] will 
be used to calculate these unknowns. The periodic disturbances are determined by calculating the 
power spectral density of the optimal control for the motion ot a spacecraft around Sun-Earth L 1 . 
An optimal control would keep the spacecraft right on the reference orbit using the least amount ot 
fuel. [Hoff93] substitutes x'and x' with the reference states to calculating the optimal u 7 . e.g. 
solving 

x" J =Tx"'+B'u'+a l4rci +a / , (3-7) 

and for B/vl 

B J u 7 =x Rj -A’x Rl -a^ -a,* (3-8) 

0xi only contains the nonlinear terms that differentiate the simulation from the reference orbit. The 
power spectral density of fl'u' then gives the values for the frequencies of these nonlinear terms. 
Here a SO i pre and a FB are the accelerations of solar radiation pressure and third body. 

3.3 - Autonomous Navigation 

The control law described above requires full state feedback, i.e. the control is a linear function of 
the state vector. In practice, the state vector is not available, and instead all one has is an estimate 
of its value generated by an observer. The Kalman filter is the optimal observer, and when coupled 
with the LQR controller, the resulting feedback system is called a linear quadratic Gaussian, or 
LQG controller. 

Autonomous navigation initiatives are underway at GSFC that cover spacecraft orbits that are 
beyond the regime in which use of GPS or the Tracking and Data Relay Satellite System (TDRSS) 
is feasible. GSFC is assessing the feasibility of using standard satellite attitude sensors, 
communication components, Doppler, crosslink systems, and GPS transceivers to provide 
autonomous navigation for missions that encompass libration points, gravity assist, high-Earth, and 
interplanetary orbits. The absolute and relative navigation solutions computed by the onboard 
Kalman filters can then be used in the LQG control of each spacecraft in a formation. 

3.4 - Decentralized LQG Design 

The model may also be cast into the form of a discrete-time distributed system, in which each 
satellite is assumed to locally process only its own measurements and generate its own controls, at 
discrete intervals. Representing the discrete time interval by the subscript i , and the sampling time 

by At, the re-cast model is: 

x, = <d,x m + £{a' + w } (3-9) 

/= i 

yj = H/x, + v/ (3-10) 

Here, x, <f>. and w are the state, state transition matrix, and process noise for all of the spacecraft in 
the formation, and 

-A! 

O, = 0(f, + At,t,) = e A&l , A'u' = JflV(0dr (3.1 1) 

t, 

The matrix B J is the same as the matrix ZTused previously, but augmented with zeros in the 

appropriate locations so as to make it compatible with the state vector containing the states of all the 
spacecraft in the formation. The process noise is used in the Kalman filter design to accommodate 
model uncertainties, and y,' and v/ represent the measurements and measurement noise of the 
autonomous navigation system at each satellite. The covariance of the process noise is W and the 
covariance of the measurement noise is V s . 

A centralized solution to estimating the state for the system above would consist ot the Kalman 
filter. The globally optimal state estimate x, for the centralized Kalman filter depends on the 




measurements from all of the spacecraft in the network. Computing x at each spacecraft therefore 
requires that all the spacecraft have all the measurement information (from all the other spacecraft), 
which in turn requires the transmission of large quantities of measurement data. ^ In contrast, 
Speyer’s approach minimizes the requirements for data transmission. Each spacecraft computes a 
globally optimal state estimate using only local measurements, along with information states an 
data vectors that compress the past measurement history into two vectors with the dimension of the 
state The approach is based upon a decomposition of the state such that: 

X = x D + x, c ( 3 - 12 ) 

where x, D depends only on the measurements (data), and x, c depends only on the controls. Based on 
this decomposition, the local Kalman filters are modified to process only that part of the 
measurement that depends on the data-dependent part of the state. 

y/=y;-H/xf ( 3 ‘ 13 ) 


The control-dependent partition propagates according to: 

x^ = O xjl, + ^ A’ u'_, (3-14) 

7*1 

with the initial condition: xf = x, . Then the data-dependent partition can be propagated and 
updated with all the data using the state transition matrix: 

x° = O x“, ( 315 ) 

with initial condition: x“ = 0 . The data dependent part of the state estimate uses only locally 

available information at each spacecraft: 

X D ' = X D ' + K'[y/ - H'x, Dj ] ( 316 ) 


where k; is the locally optimal filter gain 

K/ = (p/)(H/) t (V/)"' (3T7) 

computed using local covariance P/ , based only on local data. Note that x? is globally sub- 
optimal but locally optimal. [Spey79] shows that the globally optimal state can be reconstructed 
locally, if each spacecraft transmits its information state, (p/) i* , and a data vector, h' , to all the 
other spacecraft. Then, the globally optimal state estimate is given by. 

X° = £[(p,XP/)> + hf] ( 31g ) 

The data vector, which has the same dimension as the state, is maintained at each spacecraft using. 

h; = f h' , + g; x 1 ” ( 3 - 19 ) 


where: 
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Note that both local and global covariance analyses must be performed to generate P, and P/ . In 
some applications (e.g. those with linear time-invariant models or those with a fixed a prion 
reference trajectory), the covariance analysis can be performed offline. 


The globally optimal LQR control can be computed at each spacecraft by exchanging a set of 
vectors a'' , given by 

a" =(A') f 5,„o[(P,)(P/)'i D ' + h'l 


(3.22) 


which have the same dimension as the control, among the spacecraft. Each spacecraft transmits a 
to and receives of' from the other spacecraft, then computes the control according to 
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Note that the information exchange need only occur when a maneuver is to be performed. It the FQ 
algorithm (described below) is used, then instead the vectors h ' , which have the dimension of the 
state vector, must be exchanged at every measurement update epoch, and the globally optimal state 
is reconstructed. 


3,5 - GSFC Controller Application 

The GSFC Formation Flying Algorithm, FQ [Folt98] is adaptable to generic formation flying 
problems and permits foil closed-loop three axis orbital maneuver autonomy onboard any 
spacecraft. This algorithm, to be demonstrated on the Earth Observer- 1 formation flying mission, 
solves the position maintenance problem by combining the boundary value problem, initial and 
target states, and Battin’s ‘C*’ matrix formulation to construct a state transition matrix. In this 
example, the goal of the algorithm is for a spacecraft to perform maneuvers which cause it to move 
along a specific transfer orbit. The transfer orbit is established by determining a path which will 
carry the spacecraft from some initial state, (r 0 , v 0 ), at a given time, t 0 , to a target state, (r t , v t ), at a 
later time, ft. The target state found will place the spacecraft in a location relative to the control 
spacecraft so as to maintain the desired formation. Back propagating the target state to find the 
initial state the spacecraft would need at time t 0 to achieve the target state at time t t without 
executing a maneuver gives rise to the desired state, (rj. v<j) at time t 0 . The initial state can now be 


differenced from the desired state to find: 
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The original application of the FQ algorithm used a state transition matrix calculated using 
universal variables and the F&G series in a two-body formulation. For the application here, we 
derive the state transition matrix using the matrix exponential as shown above. That state transition 
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Where the starred quantities are based upon the position/velocity partitions of O(t 0 ,ti), and unstarred 
quantities are based on a <D(ft,t 0 ), which Battin calls the guidance matrix and navigation matrix 
respectively. If a reversible Keplerian path is assumed between the two states, one should expect 
the forward projection of the state from t 0 to ft to be related to the backward projection of the state 
from ft to to. From these sub-matrices, a C* matrix is computed as follows: 

cT »)= v -(< 0 ) k(Or 0-26) 


The expression for the impulsive maneuver applied herein follows immediately: 

AV = [c’(O]5r 0 -5v 0 


(3.27) 


4 - RESULTS AND DISCUSSION 

The centralized and decentralized LQ controllers were used with the simulation described above, 
for a two spacecraft formation. The results are shown in Figures 3 and 4. Figure 3 depicts the 
control effort in terms of velocity increments, AV. The AVs displayed are for both spacecraft which 
follow their own reference orbit so that formation flying requirement can be met. While both 
spacecraft initially started on their respective orbits, the formation requirement is to maintain a 
separation of 1km. The initial offset between the spacecraft is also approximately 1km. The peak 
shown during the first 200 days is a result of the Q and R matrices chosen for the LQR/LQG while 
the periodic oscillations indicate that there are non-linear components and perturbations not 
accounted for in the disturbance accommodation. Continuous AVs expended during the first 200 
days are less than lmm/s in each velocity component. 
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Fig. 3: AV to maintain libration orbit around Sun-Earth LI without perturbations 

Figure 4 shows the motion of spacecraft 2 relative to spacecraft 1, in a local coordinate system 
centered on spacecraft 1 . This coordinate system is similar to a radial, in-track, cross-track system, 
in that it has its x-axis along the ray connecting L\ to spacecraft 1, its z-axis along the direction 
given by the cross-product of the x-axis and the velocity of spacecraft 1 relative to LI, and the y- 
axis completing the triad. One can see that the control requirement of 1km was maintained during 
the five-year simulation. 



Fig. 4 : Motion of spacecraft-2 with respect to spacecraft- 1, in local coordinates 

To demonstrate the versatility of the GSFC formation flying algorithm, control of libration orbits 
using initial states generated by the above state space model is investigated. The orbit selected is 
shown in Figure 2. The input to the GSFC algorithm is the partitioned segments of O, the initial, 
target, and desired states as described in section 3, and the time from the maneuver to the encounter 
of the target state. The initial state is displaced in the LI X direction and LI y velocity. The only 
maneuver constraint used is that the maneuver should occur immediately at the current position. 
This placed the maneuver on the X-axis of the rotating LI coordinates. To make the targeting more 
complex, the target state also includes a position at various temporal displacements to show the 
effects on all velocity components. The position data passed into the GSFC algorithm can also be 
the information states or data from the decentralized method. Table 1 shows the AV and 










information for each of the cases studied. The time interval represents the coast time from 
maneuver to achieving the target position with the AV components in the LI X-axis and Y-axis 
directions. A seen in the table, the immediate x-axis AV required to achieve the target goals is 
proportional to the x state displacement while the y-axis AV is basically equivalent to the initial 
velocity displacement. 


Table 1 - AV (m/s) to Achieve Target State in LI Rotating Coordinate Frame 


State Offset 


10 Sec 

60 Sec 

1 Day 

30 Day 

X = 1km, y = lm/'s 

x (m/'s) 

-100 

-16 

1.2e-2 

-5.7e-4 


y (m/s) 

0.9 

0.9 

0.9 

0.9 

X = 10km, y = 1 m/s 

x (m/s) 

-1000 

-17 

1.2e-l 

-5.6e-3 


y (m/s) 

1 

1 

1 

1 

X = 10km, y =10m/s 

x (m/s) 

-1000 

-17 

1.2e-l 

-6.0e-3 


y (m/s) 

10 

10 

10 

10 


5 - CONCLUSIONS 

This work demonstrates the feasibility of using decentralized control methods to maintain a 
formation of satellites at the Earth-Sun co-linear libration point. Integration of decentralized 
control concepts within the GSFC GNCC’s formation flying algorithm is a feasible technology for 
autonomous orbit determination and control. As science and mission requirements for virtual 
satellites and formation flying missions become more demanding, the benefit/cost trade of a 
decentralized architecture versus traditional approaches becomes promising. In order to advance 
this technology to the state needed to support such missions, it is essential that this innovative 
technology could be demonstrated on future missions. The application of this LQR, decentralized 
control (LQG), and the FQ maneuvering algorithms is unlimited and can be used to fully explore 
the NASA missions outside the Earth-Moon environment. 
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